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H.T. Huynh 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

We study the numerical solutions of ordinary differential equations by one-step methods where the 
solution at t n is known and that at t n+ \ is to be calculated. The approaches employed are collocation, 
continuous Galerkin (CG) and discontinuous Galerkin (DG). Relations among these three approaches are 
established. A quadrature formula using s evaluation points is employed for the Galerkin formulations. 

We show that with such a quadrature, the CG method is identical to the collocation method using 
quadrature points as collocation points. Furthermore, if the quadrature formula is the right Radau one 
(including t„ +1 ), then the DG and CG methods also become identical, and they reduce to the Radau 11A 
collocation method. In addition, we present a generalization of DG that yields a method identical to CG 
and collocation with arbitrary collocation points. Thus, the collocation, CG, and generalized DG methods 
are equivalent, and the latter two methods can be formulated using the differential instead of integral 
equation. Finally, all schemes discussed can be cast as .v-stagc implicit Runge-Kutta methods. 

1.0 Introduction 

Collocation is an idea widely applicable to numerical analysis. In the case of numerical solutions for 
differential equations (or time-stepping schemes), for the one-step methods where the data u„ at t n is 
known and the solution u n +\ at t n +\ is to be calculated, the collocation approach can be formulated as 
follows (e.g., Hairer, Norsett, and Wanner 1987, Lambert 1991). The solution is first approximated on 
\t„, t n+ 1 ] by a polynomial P of degree 5 (for .v-stagc) interpolating the solution values at s points on 
\t n , t n \ ] ] called collocation points together with the value u n at t„. The polynomial P is determined by 
requiring that it satisfies the differential equation at the s collocation points. The solution w„ +1 is given by 
P(t n +i). For these methods, their accuracy and stability are determined by the choice of collocation points. 
For example, if the s points are chosen to be the Gauss, Radau, or Lobatto points, then the resulting 
method is accurate to order 2s, 2s - 1, or 2s - 2, respectively. Collocation methods were studied in 
(Cooper 1968, Axelsson 1969). Wright (1970) showed that the collocation process leads to an s-stage 
implicit Runge-Kutta (IRK) method. His proof will be reproduced and utilized here. 

The Galerkin method was introduced in 1915 for the elastic equilibrium of rods and thin plates (Fletcher 
1984). It was employed to solve ordinary differential equations by Hulme (1972). An introduction to both 
continuous Galerkin (CG) and discontinuous Galerkin (DG) methods for differential equations can be found 
in (Eriksson et al. 1996). The CG method seeks to approximate the solution by a continuous function which, 
on each interval \t n , t n+ j], is a polynomial of degree .v. This polynomial is determined by requiring that on 
\t n , t n+ 1 ], the weak form of the differential equation holds for all test functions that are polynomials of degree 
s vanishing at t„. Hulme (1972) showed that if an ,v -point quadrature formula is employed, then the resulting 
CG method is equivalent to a collocation method provided that the step size is bounded by certain norms to 
ensure the uniqueness of both solutions — a condition which will be removed here. (Solution uniqueness is 
not always available, e.g., for the three-dimensional Navier-Stokes equations, the problem of existence and 
uniqueness of the solution is still open.) 

The discontinuous Galerkin method (DG) is currently popular for the spatial discretization of 
conservation laws (see the review paper by Cockbum, Kamiadakis, and Shu, 2000). Formulated for 
differential equations by LeSaint and Raviart (1974), the DG method seeks to approximate the solution by 
a function, which can be discontinuous across t n , and is a polynomial of degree k on each [t„, t„+ 1 ]. At each 
t n where the solution is discontinuous, the value chosen is that just to its left — for conservation laws, such 
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a choice is called ‘upwinding’; it serves the purpose of adding numerical dissipation and results in a more 
stable method. Here, after an integration by parts, this upwind value is employed to evaluate the boundary 
term. The polynomial of degree k representing the solution is determined by requiring that on \t n , t„+ j], the 
weak form of the differential equation (after the integration by parts) holds for all test functions that are 
polynomials of degree k. Using a quadrature formula with k + 1 evaluations including an evaluation at the 
left boundary t n , LeSaint and Raviart showed that the DG formulation results in a (k + l)-stage implicit 
Runge-Kutta (IRK) method accurate to order 2k +1 or less. In addition, they proved the strong A-stability 
property (see also Bauer 1995). The method was generalized by employing the boundary values to the 
right of t n+ 1 in (Delfour, Hager, and Trochu 1981). The relation between DG and collocation methods, 
however, has not been established. On a different but related subject, it was proven in (Adjerid et al. 

2002) that for conservation laws, the DG method is superconvergent to order 2k + 1 at the “downwind” 
boundary of each cell. Concerning the basic formulation, it was shown in (Huynh 2007) that for 
conservation laws (on a quadrilateral mesh), the DG method can be formulated using the differential 
form, and the result is a simple and economical algorithm. 

In this paper, we first prove that if an 5-point quadrature formula is employed, then the CG method 
using polynomials of degree 5 is identical to the collocation method using the 5 quadrature points as 
collocation points; in other words, the condition on the step size being small enough in (Hulme 1972) is 
removed. Our proof is constructive; in addition, it shows the equivalence of the integral and differential 
forms: with appropriate choices of basis functions, one set for the space of trial solutions and another for 
the space of test functions, the CG (integral) formulation is shown to result in a collocation (differential) 
formulation. In contrast, a typical CG formulation employs (essentially) the same basis functions for the 
trial and test spaces. Next, we show that if the quadrature formula is the right Radau one (including the 
right boundary then the DG and CG methods also become identical, and they reduce to a collocation 
method called Radau HA. Compared to the proof of the fact that the DG method can be cast in the form of 
IRK by LeSaint and Raviart, our proof is more direct and leads to a specific member of the IRK class, 
namely, Radau 11A. In addition, it results in a formulation of DG using the differential instead of integral 
equation. Such a formulation can simplify the time discretization of the space-time DG scheme (for 
standard space-time DG methods, see (Van der Vegt and Van der Ven 2002)). Finally, we generalize the 
DG formulation in a manner that the resulting method becomes identical to CG. Our approach to this 
generalization does not involve the value to the right of t n+ \, therefore, it is simpler than the approach of 
Delfour, Hager, and Trochu (1981). 

Most papers on this subject are written in a highly concise manner. Often, readers can find the 
motivation and meaning of a technique or an equation only after plowing through complicated algebraic 
expressions. Such conciseness may make for an elegant style; however, it can sometimes cause 
misunderstanding. For example, in (Delfour et al. 1981), it was stated that their generalized DG method 
has the property of “superconvergence of order 2k + 1” where k is the degree of the discontinuous 
piecewise polynomial. Concerning the CG method discussed by Hulme (1972), they stated: “Note, 
however, that these continuous approximations have order 2k at the mesh points instead of 2 k +1”. It will 
be shown here that the CG method is, in fact, more accurate than DG: if s is the number of stages for the 
resulting IRK method, then, concerning highest attainable accuracy, CG is of order 2s, whereas DG, order 
25 - 1. To put it differently, for highest accuracy, CG corresponds to the Gauss quadrature, whereas DG, 
the right Radau one; as a result, CG is more accurate than DG. Note that for stiff problems, an even more 
critical criterion is stability, and here, the DG or right Radau collocation method is more advantageous. 

This paper is written in an expository manner since several different methods are involved, and a 
typical reader may be unfamiliar with one or more of them. Another goal of the expository style is to 
avoid misunderstanding. The paper is organized as follows: the collocation method is discussed in Section 
2.0; CG in Section 3.0; and DG in Section 4.0. A brief review of Legendre and Radau polynomials and a 
few examples of collocation methods can be found in the Appendices. 

We now set up the problem and introduce notations and techniques common to all methods. Note that 
the methods discussed here can be applied to systems of equations; for simplicity of notation, we deal 
only with the scalar case. Consider the ordinary differential equation (ODE) 
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u(t)=f(t,u(t)) 


( 1 . 1 ) 


with the initial condition 


u(to) -uq. (1.2) 

Let h be the step size and t„ = + nh where n = 0,1,. . ,,N. Recall that a one-step method uses one starting 

value for each step; i.e., the data u„ at time t„ is assumed to be known; the method provides a solution u n + 1 
at t n +\. For n = 0, u n is the initial condition u 0 in (1.2). Note that these methods can be applied to a variable 
step size h„; the assumption of constant step size is only for convenience. 

The following two well-known special cases are illuminating. For the first case,/ depends only on t: 

u'(t) = f(t). (1.3) 


The exact solution is 


If u n is known, the exact w„ +1 is 


u(t) = uq + [ f(x)dx . 
J t(\ 


Un+ 1, exact 


Uy, + 


[ tn+l f(t)dt. 

J tv, 


(1.4) 


(1-5) 


Thus, each one-step method results in a quantity u „ , , - which is a quadrature formula approximating 
the integral above. 

For the second case,/= A .u where A is a complex constant; therefore, 


u'(t) = A ,u(t) . 


The exact solution is again obvious: 


u(t) = UQ e K ( l 


If u„ is known, the exact u n +\ is given by 


u n+\, exact 


U n 


e Xh . 


( 1 - 6 ) 


(1-7) 


( 1 . 8 ) 


Each one -step method yields a solution u n +\. Define the stability function R by 

u„+i = u „R(kh). (1.9) 

Then with z = A h, (1.8) implies R(z) approximates e . If the method is of order p, then the local error is 

E(z)= e z -R(z) - 0(zP +1 ). (1.10) 


In other words, 


R(z) - 1 + z + — + ... + — + 0(zP +1 ) . 


2 ! 


P- 


( 1 - 11 ) 
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The converse, however, does not hold: due to nonlinear errors, (1.10) or (1.1 1) does not imply that the 
method is of order p. (Note that the quantities 0(z p ' 1 ) in (1.10) and (1.1 1) are different from each other). 
The stability domain is 


S = { complex number z such that |i?(z)| < 1 } . 

With z = Xh, the solution for (1.6) after n steps is u„ = Ui, R(z) n . If z is in S, then \R(z)"\ < 1 ; therefore, 

\u„\ < ]«o| for all /;, i.e., the solution is bounded for all time. Next, a method is A-stable if the 
corresponding S contains the left half of the complex plane: 

S =3 { z ; Re(z) < 0 } . (1.12) 

An A-stable method (such as the trapezoidal rule (B.3) below) can have the following property, which is 
not always desirable, 


lim R(z) = 1 . 

Z — > 00 

With such a property, for an exact solution that damps quickly (say, A 1000 '), the approximate solution 
damps very slowly. A more desirable property is L-stability: a method is L-stable if it is A-stable and 

lim R(z) - 0 . (1.13) 

Z — ^ GO 


L-stability implies that the method is suitable for stiff problems (the X values or eigenvalues of a stiff 
problem have magnitudes in a wide range, from very small to very large). 

The following rescaling technique is employed extensively below. Denote /„ = \t n ,t n + 1 ] and I = [0,1]. 
Instead of it is often more convenient to work with I. For t on set 

x = (t-t„)/h . (1.14a) 

Then x varies on I. The inverse maps I onto 

t = t n + x h . (1.14b) 

Each function g(t) on /„ corresponds to a function g(x) on /, namely, g(x) = g(t n + x h ) . Here, we use 
the notation g(t„ + xh) instead of g . Denoting g 1 = dgldt, we have, by the chain rule, 

~g(t n +xh) - h = h g\t) . (1.15) 

dx dt 

As for integrals, 

[ " +> g(t)dt = h f g(t n +xh)dx . (1.16) 

j t n jo 

2.0 Collocation Methods 

To describe these methods, let c„ i= 1,. . ,,s be (collocation) points in ascending order on 7, 

0 < Cj < 1 and c, < Cj for i < j . (2. 1 ) 
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Let the collocation points on I„ be defined by 


7/z, i ~ Oz Cf h . 


(2.2) 


Suppose, for the moment, the solution values u n j at t„j, i = 1,...,5, are known. These 5 values together with 
the value u„ at t n determine a polynomial P = P(t) of degree s (the case c x = 0 will be discussed later), 

P(t n ) = u n (2.3) 


and, for i = 


P(tn, i) ~ W/z, i • (2-4) 

The quantities P'(t n i ) = (dP/dt)(t n i ) and f(t n j,u nJ ) can then be evaluated for i = 1 ,. . ,,s. The collocation 
method seeks a polynomial P that satisfies the following implicit equations: for i = 1 ,. . .,s, 

P\tn, i ) = fih z, i , M„, i) = f(t n , i , P(t„, ,)) (2.5) 

Once P is determined, the solution at t„ +l is given by 

^zz+l — Pfn+l) ■ (2.6) 

Two remarks are in order. First, the approximating polynomials for u and / of the ODE are of different 
degrees. Indeed, if u is approximated by P of degree s, then u' is approximated by P' of degree s - 1. 

Since u' =f we wish to approximate /by a polynomial of degree s - 1. For the collocation method, this 
polynomial is determined by the values of f at t„ Al . ..J, LS (and not the value at t„). As an example, for the 
equation u' = Xu, the function u of the left hand side is approximated by P, whereas, u of the right hand 
side, by the values at the collocation points t nA ,. . ,,/ s (a polynomial of degree s - 1). 

The second remark concerns the case c, = 0. Flere, the derivative P'{t nA ) = P'(t n ) =/ (t n ,u n ) can be 
calculated explicitly. Therefore, P is determined by (2.5) with i = 2,...,s and 

P{t n ) = u n and P\t n ) = f(t n ,u n ). (2.7) 


2,1 Collocation and Implicit Runge-Kutta (IRK) Methods 

It was shown by Wright (1970) that the above collocation method results in an .v-stagc IRK method. 
The proof of this fact (Lambert 1991) is reproduced below since it will be employed for our proof of 
equivalence between CG and collocation methods. It amounts to expressing P' in terms of certain basis 
functions and then integrating P' to obtain P in IRK form. 

Consider the s collocation points c u . ,.,c s on [0,1]. The values (of a function) at these points determine 
a polynomial of degree 5-1. With i fixed, 1 < i < s, let L,(x) be the Lagrange polynomial of degree s - 1 
defined by L,(Cj) = 5,, for / = l,...,s; i.e., L, takes value 1 at x = c, and 0 at all other Cj,j / i (see Fig. 2.1), 


fl 

7=1, j*‘ 



(2.8) 


Then L h 1 < i < s, form a basis for the space of polynomials of degree 5 - 1 on [0,1]. 
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Figure 2.1. — (a) Functions Lj,j = 1,2,3 for the case of 3 Gauss points (s = 3) as collocation points on the interval 
/ = [0,1]. (b) The same functions for the case of 2 right Radau points (C 2 = 1). 

For i = l,...,s, set 


k i ~ J ~(tn, i 5 M/i, /) • 


(2.9) 


Next, observe that P' = dP/dt is a polynomial of degree 5 - 1 . Using / for the index instead of i (the reason 
for this switch will be clear in (2.1 1)), since P' interpolates the s data points (t„ + Cjh,kj), 


P'(t n +Th) = Y J j =l Lj(P) k 


v • 


( 2 . 10 ) 


Concerning the above left hand side, with t = t„ + xh, we have 


f C ' P\t n + x h) dx = [ tn +C ‘ h P\t) dt = P{t n + Ci h) - P(t n ) . 

JO J tv. 


Equation (2. 10) then implies, for the z-th stage, 


p(t n+Ci h)-p(t n ) - f/ (x ) 


dx 


kj. 


(2.11) 


As for the solution u n+l , 


M/i+l 


= P(t n +h)-P(t n ) = S fj 1 X y (T) 


dx 


k j ' 


( 2 . 12 ) 


Motivated by (2.1 1), set 


f c > 

J Lj (t) dx , 


(2.13) 


and, by (2.12), set 


b J = Jo 


dx . 


(2.14) 


Then, with u n j = P(t n j ) and kj =f ( t n j , u n J), (2. 1 1) implies the following IRK form of the collocation 
method: for i = l,...,s, 
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(2.15) 


^ n. i — Mn b i | @ij kj . 

After obtaining u n j (and k,) by solving the above system of s implicit equations, the solution is given by 

u „+ 1 =u„ + h^ S bjkj. (2.16) 


This completes the proof. 

The Butcher array of an IRK method consists of c„ a, p and bj arranged as follows 


C\ 

a \ i 

d\ s 

c s 

a s\ 

. . s 


b\ 

.. b s 


Here, for each i - th row, a p j = 1 ,. . ,,s are the weights of a quadrature formula as will be shown in (2. 1 8). 

Note that if c x = 0, then t nA = t m u nA = u n , and ay = 0 for all j; in other words, the first row of the 
Butcher array is identically zero. In addition, the quantity k x =f ( t„,u n ) can be calculated explicitly. Since 
u nA = u„ is known, (2.15) for i = 2 then yields the equations to calculate m„ i2 ,. . 

Concerning the stability function for the IRK method, define the s><s matrices A = [a,,] and 
I = identity matrix, column vectors of s entries b = [/u,... A] 7 and e = [1,1,. . .,l] r where the superscript T 
stands for transpose, and det = determinant. Then, for the IRK method, the stability function R can be 
calculated by one of the following two formulas (e.g., Lambert 1991) 

R(z ) = 1 + zb T (l - zA) -1 e 


or 


det [1 - zA + zeb T ] 
det [I - zA] 

2.2 Quadratures Associated with IRK Method 

The above IRK method relates to quadrature formulas as follows. By (2.16) and the definition of kj, 

u n+\ ~ u n — b ^ i bj U n j . 

The corresponding quadrature formula is, for any continuous function v on 7= [0,1], 

j o v(x)dx ~ X /=i b J v(c j) ■ ( 2 - 17 ) 

Here, bj are the weights given in (2. 14), and the collocation points Cj are the evaluation points. Similarly, 
by (2.15), 


= hY s , 

^7=1 


®ij u n,j • 
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The corresponding quadrature formula is 

jj' v(x)c/x ~ Yjj ={ a ij v ( c j)- ( 2 - 18 ) 

Here, we use the s collocation points on [0, 1 ] to evaluate the integral from 0 to c,. 

An observation concerning accuracy of these quadratures is in order. With appropriate choice of 
evaluation points on/, formula (2.17) has a degree of precision of up to 2s - 1. (Recall that a quadrature is 
of degree of precision m if it is exact for polynomials of degree m or less.) Formula (2. 1 8) for the stages, 
however, has a degree of precision no higher than s - 1 since special points (say, Gauss points) on [0, 1 ] 
are generally not special on [0,c,]. For example, if we use 2 Gauss points on [0,1] as collocation points, 
then the degree of precision for (2.17) is 3, but that for (2.18) is only 1. As will be shown, the resulting 
collocation method is of order 4. 


2.3 Basis Functions Lj 

We next introduce the basis functions Lj , which will be employed in the proof of equivalence 
between the collocation and CG methods. With t = t n + xh, similar to (2.1 1), by (2.10), 


P(t)~P(tn) = h lL] = Ml L J^ 


k J- 


(2.19) 


Denote, for / = 1 , . . . ,s. 




(2.20) 


Then L: is of degree s since Lj is of degree s - 1. Some additional noteworthy properties of Lj are: 


d L / / dx - L ; 


(2.21a) 


in addition, 


^ ■ / ( 0 ) 0, Lj(Cj) cijj , and Lj(X) bj , 


moreover, since 



(2.21b) 


Z - 


x. 


Next, by (2.19), 


(2.22) 


P(t) -P(t n ) = P(t n + xh) - P(t n ) = h^ s J=l kj Lj (x) . (2.23) 
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Figure 2.2. — (a) Functions Lj j = 1,2,3 for the case of 3 Gauss points on / = [0,1]. Note that dLj(x)/dx = Lj(x ) ; 
therefore, the graph of Lj has slope 1 at c, and slope 0 at all other c, as can be seen by the slopes at the 
dots, (b) Functions Ly for the case of 2 (right) Radau points. 


For convenience, set 


L 0 = 1 (2.24a) 

and 

A'o = P(t n )/ h = u„ / h . (2.24b) 

Then P(t n ) - h Atq Lq , and (2.23) can be expressed as 

P(t n +xh) - ^z; 0 kjLjix) . (2.25) 

Observe that Lj ,j = 0 arc .v + I polynomials of degree no higher than ,v. We will show that they 
are linearly independent; thus, they form a basis for the space of polynomials of degree s or less. That is, 
we wish to show that if 


Z 7 Lo a J = 0 ’ (2.26) 

then 

a y =0 for j = 0,...,s. (2.27) 

To this end, observe that (2.26) implies Z --o 01 j (®) • Since Lj( 0) = 0 for j = l,...,x, we have 

ap Tq( 0) = 0 . By definition, Lq = 1 ; as a result, 


ao = 0. 


Equation (2.26) then takes the form 



Note the starting value of 1 for j. Differentiating the above, we have 


(2.28) 


(2.29) 
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Xy=i a j L J=°- 

Since Lj,j = 1,. . .,5, are independent, we conclude that a.\ = ... = a s = 0. This fact and (2.28) complete the 
proof that Lj are independent. 

The above observation implies that Lj ,j= 0,. . ,,s, form a basis for the space of polynomials of degree 
no higher than s. Thus, if Q is any polynomial of degree s or less, then Q can be expressed as a linear 
combination of Lj\ 


Q(t n +xh) = hYjj=o a J M x ) (2.30a) 

where the coefficients a , relate to Q by 

«o =Q{t n )lh and aj = Q\t nJ ) for j = . (2.30b) 

Also note that if Q is a polynomial of degree s and Q(t„) and Q'(t n J),j = 1 ,. . .,.s\ arc known, then Q is 
given by (2.30). Finally, examples of the Gauss, Radau, and Lobatto collocation methods can be found in 
appendix A. 

3.0 Continuous Galerkin (CG) Methods 

The CG method seeks to approximate the solution by a continuous function which, on each interval 
\t,„ t„+ 1 ], is a polynomial U of degree s determined by using the weak form of the differential equation. 
Again, assuming that U(t„) = u„ is known, we wish to calculate u n + 1 = U(t n + 1 ). 

We need some preparations. Let [a,p] be any interval; here, it is either 1= [0,1] or /„ = [t n , t n +{\. For 
simplicity, unless otherwise stated, we deal only with continuous functions (in fact, polynomials) on 
[a,p]. The inner product of two functions Vi and V 2 is defined by 

(vi,v 2 ) = \^vi{t)v 2 {t)dt . (3.1) 

J a 

Next, denote by P s [a,p] the space of polynomials of degree 5 or less on [a,p]. In addition, denote by 
P 0 ' [a,P] the subspace of P s [a,p] consisting of polynomials that vanish at the left boundary a. When 
there is no confusion concerning the interval, we use the notation P and P ( ) . Note that P' is of dimension 
s + 1 , and Pq , dimension s. 

In general, the integrals, e.g., the inner product (3.1), are carried out by approximations rather than by 
exact integration. We can employ the quadrature (2. 1 7) on / (the Cj are now the evaluation points): 

Jo v (t></t * Z -= iV^y)- ( 3 - 2 ) 

The quadrature on /„ differs from the above by a factor of h: for any continuous function v(t), 

f " +l v(t)dt * h ( 3 ‘ 3 ) 

~ tn J 

We often use (3.2) instead of the above (i.e., the factor h is understood) when there is no confusion. 

The trial space or space of trial solutions on /„ is defined by 

S = {polynomial U of degree s such that U (t n ) = u n } . (3.4) 
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In spite of its name, S is in fact not a space, but is a hyperplane in P \t n , t n+] \. 

A test space or a space of test functions on /„, commonly denoted by V, is a subspace of P which has 
dimension .v. Among the most commonly used test spaces is the space consisting of polynomials that 
satisfy the homogeneous boundary condition at t„ 

V = {v of degree s such that v(t n ) = 0} , (3.5) 

i.e., V = Pq . Note that if JJ is in S, then U - u n is an element of the above V. 

Next, the weak form of the equation u' =/ can be written formally as 

(«',v)= (/». (3.6) 

The CG method seeks a solution U in the trial space S that satisfies, for all v in the test space V, 

(U',v) = (/» . (3.7) 

In other words, the projection of U and / ( t,U(t )) onto V are the same. 

We now show that the test space, as opposed to the case of trial space, must be a subspace, and it must 
have dimension s. Indeed, first, suppose (3.7) holds for functions V\ and v 2 (which take the place of v). Let 
a and P be real numbers. Then, one can easily verify that (U\ aiq + Pv 2 ) = (f avi + Pv 2 ); thus, the test 
space must form a subspace. Next, since U is of degree s, and one condition is known, namely, U(f„) = u,„ 
s conditions remain to be determined. Therefore, the test space must be of dimension ,v. 

We next discuss the choice of test spaces. Instead of solving for U, we can solve for U— u„. Since 
U - u„ is in V, and V defined by (3.5) has the correct dimension, it is sensible to use V as a test space. 

The following argument results in another choice of test space. Since U is of degree s — 1 , for the two 
sides of (3.7) to match, we should approximate /by a polynomial also of degree 5-1. Such an 
approximation can be obtained by the projection of / on P s_1 . For this projection to be identical to U, it 
suffices to require that (3.7) holds for all v in P , i.e., the test space be P' 1 . It was observed in (Fletcher 
1984) and also in (Eriksson et al. 1996) that this test space yields a more accurate solution than the test 
space Pq . We will show in (3.34) below that with the quadrature (3.2), or equivalently, (3.3), the two test 
spaces V = Pq and V = P s_1 yield identical results. The following conclusion can then be drawn: for the 
CG formulation with the test space V = Pq , the method using Gauss quadrature is more accurate than that 
using exact integration. (This fact is contrary to the common belief that more accurate integration 
formulas yield more accurate solutions.) 



Figure 3.1 . — (a) Basis functions <j>y j = 0,...,3 for the case of 3 Gauss points on / = [0,1]. (b) Basis functions 4>o, 4>i , and 
§2 for the case of 2 right Radau points. Note that <]>s form a basis for Pq . 


N ASA/TM— 20 11-216340 


11 


In addition to Lj defined by (2.20), we will need the following basis for P s . On /, assuming c\ > 0, set 


c 0 = 0. 


Then, c 0 > C\ > .... c s < 1. Let <t) 0 ,..., <K be the corresponding basis functions: for each j. 


<m t ) = n 

i= 0, i^t j 


X — Ci 


C J~ C ‘ 


(3.8) 


(3.9) 


That is, 4> y(c,-) = 8,, as shown in Figure 3.1. (Note that the definition of L, in (2.8) is different from the 
current definition in that it does not include c 0 ). 

Since <()o(0) = 1, we can employ ())o to deal with the boundary condition at t n . For all j > 1, <jiy(0) = 0; as a 
consequence, (|>i,. . ., <t>^ form a basis for Pq . Also note that in the case of s Gauss points, ())o defined above 
is identical to the Legendre polynomial of degree .v on /; in the case of right Radau points, it is identical to 
the Radau polynomial (see Appendix). In the case of Lobatto points, however, (j> 0 is not the Lobatto 
polynomial; this case corresponds to c x = 0 and will be discussed later. 


3,1 CG and Collocation Methods 

We claim that if the quadrature formula (3.2) is employed to evaluate the inner product, then the 
resulting CG method is identical to the collocation method with quadrature points as collocation points. 

To prove the above claim, we will start with the CG solution U and show that it can be expressed in 
collocation form. The key ingredient of the proof is the choice of {Lj}j =0 defined by (2.20) as a basis for 
S and {(()/' } J = [ defined by (3.9) as a basis for V = Pq . Since U - u„ is in V, and {Lj} j =1 is a basis for V, 

U — u n can be expressed in terms of these basis functions: 

U(t n +xh) - u n = h ^ S J=l djLj (x) . (3.10) 

As discussed in (2.30), dj and U are related by, for j = 1,. . ,,s, 

d j =U'(t n j ) . (3.11) 

We wish to show, by using the weak form (3.7), that dj =f (t n j,U(t nj )). To this end, because d Lj / dx = Lj , 
by differentiating (3.10) and employing the chain rule, 

U\t n +xh) = X - = i d J L J^) ■ ( 3 - 12 ) 

Next, recall that U is a CG solution. Noting that (j)i,. . .,(j), form a basis for V, they can serve as test 
functions: for i= l,...,s, 


(tf',4>,-)=(/,4>.-). (3-13) 

We will show that under the quadrature rule, the left hand side above yields U'(t nJ ), and the right hand 
side, /(!„,;, U(t n j)). Indeed, by (3.12), 


(U'A) = X-=i d j^ L j^i)- 


(3-14) 
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When the inner product is evaluated by quadrature (3.2), we use the notation e.g., 

(v, w) q = X/=t bi 

In (3.14), Lj vanishes at all collocation points except at c h and (j>, vanishes at all collocation points except 
at Cj. Consequently, the only nonzero term for the sum is d,{L h ty,). Since L,{c,) = <|>,(c,) = 1, employing the 
quadrature rule, (3.14) implies 


{U'A) q = b x d x . (3.15) 

Concerning the right hand side of (3.13), under the quadrature rule, only the values of / at the quadrature 
points are needed. Since /((), takes on the value / (?„,„ U( ?,,,)) at x = c, and the value 0 at all other c„ 

(fAi)q = bif(t n; i,U(t nJ )). (3.16) 

Again, under the quadrature rule, by (3.15) and (3.16), equation (3.13) implies 

d i — f ifn,i riUifn,i )) • 


From the above and (3 . 1 1 ), for i = 1 , ... ,s, 


U\t nJ ) =f(t nj ,U(t nj )) 1? 

That is, U plays the role of P in the definition of the collocation method (2.23). This completes the proof. 

3.2 Standard CG Method (Via Basis Functions) 

In (3.7) above, the CG method is formulated via trial and test spaces. It can also be formulated via 
basis functions (Hulme 1972, Delfour et al. 1981). This formulation is presented here with more 
explanations and will be employed later in (3.35) to (3.45). Let (p 0 , . . . ,(p, s be a set of basis functions for 
P S [0,1] with the following properties. First, 


<Po(0) = l 

so that cp 0 can deal with the boundary condition at t„. Next, for j = 1 ,. . . ,,s, 


(3.18) 


9 j (0) = 0 (3.19) 

so that (pi,. . ,,cp s form a basis for Pq . Note that Lj defined by (2.20) and (2.24a) and (j), by (3.9) possess 
these properties. Set 


d 0 = u n / h . (3.20) 

Then the trial solution U of (3.4) can be written as 

U(t n +xh) = h ^ =o c/ y <P/(x) (3.21) 

where d 0 is given by (3.20) and cl\,. . . ,d s remain to be determined. Taking d/dx of the above, 
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(3.22) 


U\t n + xh ) = y n dj d(pj(T)/(h . 

Using the test function v = cp, where i = 1,. . .,5, defined by (3.9) (these cp, form a basis for V), the weak 
form (3.7) implies 


Z; o (CP, d(pj / dx) dj = (<p ( 3 . 23 ) 

We thus obtain .v implicit equations (since / also depends on d,) for .v unknowns d u ... ,d s . 

Note that in general, the expression f{t, y^. n d j (p/(t)J > which depends on the unknowns dj, can be 

complicated; as a result, evaluating (cp,,/) via exact integration may become difficult. It is often more 
convenient to evaluate the inner product by a quadrature formula employing only the values / (t nJ , U(t nj )) 
at the quadrature points. 


3.3 Equivalence of Collocation and CG Methods Via Approximating the Dirac Delta 
Function 

If we use the Dirac delta function, the proof of this equivalence is shortened considerably. The basic 
idea is that on P with an inner product by quadrature (3.2), the Dirac delta function 8, has the same effect 
as c ybj where ()), is the basis function given by (3.9) and b, is the quadrature weight. Indeed, on I = [0,1], 
let 8, be the Dirac delta function at c, defined by, for any v in P , 


(v,8,) = v(c,). (3.24) 

Using the quadrature rule (3.2), again for any v in P , 

0, <t>, lbi) q = v(Cj) = (v,8,) . (3.25) 

That is, concerning the inner product on P via the quadrature rule, we have 

<th / b t - 5j . (3.26) 

The collocation method (2.5), with P defined by (2.3) and (2.4), can be written as 

(P',8,) = (/,8,). (3.27) 

Therefore, by (3.25), for / = 1,. . ,,s, 

( P'M) q = (/,< l>;V (3-28) 

The above is the CG form (3.13) with the inner product evaluated by the quadrature (3.2) and U replaced 
byP. 

Note that for the above argument to hold, <(>,■ is required to be of degree s and thus, is defined by ,v + 1 
conditions, but only the .v values of ()), at the collocation points are needed in the above proof. The extra 
condition can be arbitrary (i.e., the condition <(),(0) = 0 is not required) as discussed further in (3.33) 
below. 
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3.4 The Case C\ = 0 (and cr > 0) 


Examples for this case are the Lobatto and the left Radau quadrature points. Via the collocation 
approach, as discussed in (2.7), for this case, in addition to the boundary value u n , the derivative 
u ' n = / (x„,u n ) is also known (easily evaluated). Concerning the CG approach, it needs to be modified to 
incorporate the condition that u' n is known. To this end, the trial space is defined as 

{polynomial U of degree s such that U(t n ) = u n and U'„ = f(x n , u n ) j- . (3.29) 

This trial space is of dimension s — 1. The test space can be modified accordingly: 

{v of degree s such that v(t„) - 0 and v'(t n ) = 0} . (3.30) 

This space is also of dimension s — 1 . For the case c\ = 0, the CG method seeks a solution U in the trial 
space (3.29) that satisfies, for all v in the test space (3.30), 

(U',v) = (f,v) . 

Next, the definition of Lj in (2.20) remains valid. In addition, the definition of the basis functions ()>, in 

(3.9), except for <()o and <(>!, also remains valid. Since c 0 = c\ = 0, to modify <|)o, we define it to be a 
polynomial of degree s such that 

4>0 (0) = 1, <t>(/(0) = 0 , and <t>o ( c i ) = 0 for i = 2,...,s . (3.31a) 


That is 


where 


<t>o(x) = 0-ax) Yl 


C; -x 


i= 2 


a = 


n 


C; -X 


Ki = 2 Cf J 


( 0 ). 


As for <)>i, it is defined to be a polynomial of degree s such that 

(() i (0) = 0, (|i i'(0) = 1 , and <() i (c, ) = 0 for i = 2 . 


(3.31b) 


(3.32a) 


That is 


<t>l(x) = X fl (3.32b) 

7=2 c ' 

Note that 4> 0 serves the purpose of matching the value u n , and 4*!, the slope u n ' = f(t n , u„) . 

The claim of equivalence between CG and collocation methods still holds when c\ = 0. The proof, 
which employs the basis functions <t>2,- - . ,4)^- for the test space (3.30), and the basis functions Lj for the 
trial space (3.29), is essentially the same as that of the case c x > 0 and is omitted. 
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(a) (b) 

Figure 3.2. — (a) Basis functions <t> 3 for the case of 3 Lobatto points, (b) Basis functions <]>o,<|>i, and (j) 2 for the 
case of 2 left Radau points. Note that § 2 ,—, <t> s form a basis for the test space defined in (3.30). 


3.5 Test Spaces 

We claim that using quadrature (3.2) with evaluation points c,’s, the test spaces P 1 and Pq (or V) 
yield identical results. The reason is that on P s with such a quadrature, as far as the inner product is 
concerned, (j>, defined by (3.9) has the same effect as I, defined by (2.8) for i = 1,. . ,,s. Indeed, for any 
continuous function v, with i varies between 1 and s, 

(v ,§i)q = (v,Li)q = (v, bi&i) q = bi v(Ci ) . (3.33) 

Thus, if U is a solution of the standard CG method, then since (j), is in the test space Pq , (C/,<j>,) ? = (f,(j>,) ? . 
As shown above, we can replace by L h 


(U',L i ) q =(f,L i ) q . (3.34) 

Noting that L h i = 1 ,. . ,,s form a basis for P i_1 , the claim follows. 

Readers who are interested only in the main results can omit the rest of this section with no loss of 
continuity. 

3.6 A Dilemma 

The question then is: What causes the difference between the established fact that the test space P‘ s_1 
yields a more accurate solution than the test space Pq and the above claim that under the quadrature rule, 
the test spaces P s 1 and Pq yield identical results? To answer this question, we need to derive the 
solutions using exact integration. This derivation is similar to that in (Fletcher 1984); the key difference, 
however, is that instead of u' = u, we use the equation 


u' = Xu. (3.35) 

As a consequence, we can study not only which test space yields a more accurate numerical solution, but 
also the stability as well as order of accuracy of the corresponding method. We will show that with Pq 5 as 
the test space, an appropriate quadrature formula results in a method more accurate than that by exact 
integration. In other words, such a quadrature provides the cancellation needed for high accuracy 
whereas, with exact integration, the cancellation no longer occurs. This observation contrasts the common 
belief that a more accurate integration procedure results in a solution with better accuracy. 

To solve the above ODE, consider the following basis functions for P S [0,1] 

cp/ = W, j = 0,...,s . (3.36) 
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Clearly, cp/O) = 0 for j = 1 , . . . ,s, and cpi,..., cp* form a basis for P ( j = V . Next, set 

d 0 = u n /h (3.37a) 

and 

U=h H% o d jVjW ( 3 - 37b ) 

as in (3.20) and (3.21) where, d\,. ,.,d s are to be determined. For the equation u' = Xu, the standard CG 
formulation, namely (3.23), implies 


to, d(pj/dx)dj = ^X/=o (Vi’Vj)dj • (3-38) 

Here, for the test space Pq , i varies from 1 to s; for the test space P , from 0 to s - 1 . Moving the two 
terms corresponding toy = 0 to the right hand side and the rest to the left, we have 

X/=i <#//<*) - M<P;> 9/)] dj = — (<p,- , rfcpo / rfx) + X, (cp,- , cpo) do. (3.39) 

Since cp 0 = 1 , the first term on the right hand side above drops out: 


d<$jldi)-X(<Vi, cp j)]dj = X (cp,- , cp 0 ) • 

With exact integration, we obtain, for all i and j, 

(cp,- , c/<p ; / dx) - [ jx' + J~ l dx - 

Jo i + j 


and 


ft . . 1 

(cp/ , cp/) = L ^ l+J dx = - — — . 

J JO H- 7 + 1 


(3.40) 


(3.41a) 


(3.41b) 


Consider now the case s = 2. Assuming that u„ = 1, h = 1, then by (3.37a), c/ 0 = 1 . With exact 
integration, the two cases concerning different test spaces are as follows. 

For the test space P‘ s_1 = P 1 , using (3.41), the system (3.40) with i = 0 and i = 1 takes the form 


1_A 

V 2 3 


1-% A 


A 

47 


fdA 

V d 2 


X 

& 




J 


(3.42) 


The solutions are d\ = [6X(2 - X)]/(12 -6X + X ' ’) and d 2 = 6X7(12 - 6X + X 2 ). Since u n + \ = 1 + d\ + d 2 , 


u n + 1 


12 + 6X + X 2 . 1 + z/2 + z 2 /12 

— or R(z) = . 

12-6X + X 2 1-z/2 + z 2 /12 


(3.43) 


Thus, the method yields the same solution as the collocation method (2.37) using two Gauss points. It is 
of order 4, and the error can be found in (2.38). The stability region is shown in Figure 2.3(b). 
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For the test space Pq , using (3.41), the system (3.40) with i = 1 and i = 2 takes the form 


( 1 X 

2 X\ 

(d^ 


f 1 

2 3 

3 4 


2 

1 X 

2 X 



1 

V 3 4 

4 5 ) 


U 


(3.44) 


The solutions are d\ = [4 A, (5 - 3A)] / (20 - 12 X + 3X 2 ) , di = 10A, 2 / (20 - 12A, + 3X 2 ), and 


The error is 


20 + 8A + A 2 
Un+l ~ 20-12A + 3A, 2 


or 


R(z) 


20 + 8 z + z 2 
20-12z + 3z 2 ‘ 


E(z) - e z -R{z) = 0.016667 z 3 + 0(z 4 ) . 


(3.45) 


(3.46) 


Thus, the method is of order only 2. As |z|— xx>, R(z )—+ 1 /3; as a result, the method is not L-stable. It is A- 
stable, however. 

Note that if the Gauss quadrature is employed to evaluate the left hand side of (3.38), the result is the 
same as that by exact integration since ((p/,d(p/</T) is of degree 2s - 1 or less. For the right hand side of 
(3.38), however, when i = j = s, cp„<p / = r 1 ', and the Gauss quadrature yields a result different from that by 
exact integration. Thus, for the test space Pq , since the Gauss quadrature for CG results in a method 
equivalent to the Gauss point collocation, such a quadrature yields a CG method of order 4. For the same 
test space, i.e., Pq , exact integration, as shown above, results in a CG method of order only 2. Thus, 

concerning the CG formulation, a well-matched quadrature can provide the cancellation needed for high- 
order accuracy, whereas exact integration may not, and the result is a less accurate solution. 


4.0 Discontinuous Galerkin (DG) Methods 

The DG method seeks to approximate the solution by a function which is allowed to be discontinuous 
across each t„ and, on each interval is a polynomial of degree .v - 1 denoted by u h (i.e., one degree 
lower than that of U, the corresponding CG solution). This polynomial is determined by using the weak 
form of the ODE together with an integration by parts. Here, we assume that u h (t,~) , which plays the 
role of is known. For the first interval, «/,(t 0 ) = uq • We wish to calculate Uh{t n ~ x ) , which plays the 
role of u n +\. 

On /„, formally (the correct expression is (4.4) below), the weak form for the equation u' =/ is, for any 

v in P‘ s_1 , 


[ " +l Uh' (t)v(t) dt = [ ' +> f(t,u hit)) v(t) dt . 

" t n " 


(4.1) 


The above equation contains no effect from the initial condition u /, (t„ ) = u n . To involve this condition, 
we use integration by parts: 

f " +1 Uh (t)v(t) dt - [u h v\ Y 1 ~ [ " +1 Uh(t)v'(t) dt . (4.2) 

■’ C " ■’ tn 

At each t n , n = (),..., /V, the value u h is not well defined since u /, (t n ~ ) is generally different from u h (t„ + ) . 
By assuming that ‘time marches forward’, the value it/, (t„ ) is employed for the boundaiy evaluations. 
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Such a choice for the case of an advection (where ‘the wind blows from the left’) is called the ‘upwind’ 
choice, a term employed here as well. The choice adds numerical dissipation to the method whereas the 
‘centered’ choice, i.e., \u/, (t„) + Uh(t„ + )\/ 2 generally results in no numerical dissipation. Thus, between 
u h ( t n ) and Uh (tn ) , the upwind value is u/, (t n ~) . Note that for conservation laws, we deal with the 
upwind flux; here, the flux is u itself; therefore, we deal with the upwind value. As for the test function v 
on /„, the values v(t n ) and v(t„ i ) involve no ambiguity. The above equation then takes the form 

u h (t)v(t) dt = \u h (t n+l )v(t n+l )-u h (t n )v{t n )] - u h (t)v'(t) dt . (4.3) 

J tn J t n 


The DG formulation is the following: with the value «/, (t n ) given, we wish to find, on /,„ a polynomial 
u h in P that satisfies, for all v in P ', 

t p t J_) 

\uh(tn + \)v(tn + \)- u h(tn )v(t n )] ~ u h (t)v'(t) dt = f{t,u h {t))v{t) dt . (4.4) 

Jt n Jt n 

As in (LeSaint and Raviart 1974), we need another integration by parts. On since u h is a polynomial, 
the values at the two boundaries are respectively identical to m/,(G + ) and u/ I (t n ^ l ) . Consequently, 

u h (t)v\t)dt = [u h (t~ +l )v(t n+l )-u h (t n )v{t n )\ - u h (t)v(t)dt . 

j t„ j t„ 

Note that this equation involves no upwinding. Substitute it into (4.4), i.e., after integrating by parts twice, 
we obtain the following DG formulation: with u h (t ~) given, find «/, in P s 1 that satisfies, for all v in P v ', 

I (* ^ f \- pi # I* ^ n-\-\ 

~ \ Uh(t„ ) Ufi (t n )]v(G) + u h (t)v(t) dt - f(t,u h (t )) v(0 dt . (4.5) 

J t„ J t„ 

Between the above and the formal expression (4.1), the difference is the ‘correction’ term 
- \u h {t„) - u/, (Jn )] v 0n ) ■ This term serves the purpose of enforcing the upwind value «/, (t„ ) at the left 
boundary of The upwind value «/, (t„+| ) at the right boundary is built in as shown in (4.4). It is also the 
solution we wish to calculate. 

As an example, we solve u'(t) = 2 1 with initial condition u{ 0) = 0 via the above DG formulation. The 
exact solution is obvious: u{t) = t 2 . We use a linear element u h , i.e., 5 = 2. For simplicity of notation, set 
a = t n and b = t n+i . Assuming that u h {a~) is exact, i.e., u h {a~) = a 2 , we wish to calculate u h on the interval 
/„ = [a,h]. First, since u h is linear, the problem reduces to calculating «/,(/T)and u h {ct). Next, noting that 
h = b - a, on the interval (a, b), u h '= [ Uh(b ~ ) - ui,(a + )\/h. As a result, (4.5) implies 

1 b (• b 

- [a 2 - ui 1 (a + )\v(a) H — \_ u h{b~) -M/,(a + )] v(t) dt - 2 1 v(t) dt . (4.6) 

h J a J a 

Set v = 1 , the above results in 

- [a 2 -u/,(a + )] + [iih(b~)-Uh{a + )] = b 2 -a 2 . 

Thus, the solution at t n+ \ =b is 


u h {b ) = b 2 . 

This solution is exact. Now, setting v = t for (4.6), it follows that 
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1 


- [ a 2 -w/,(a + )] a + 


— [b 2 -u/ 7 (a + )] — ( b 2 - a 2 ) = — [b 2 -a 3 ] . 

(b-a) 2 3 


Consequently, 


«/,(a + ) 


(2a 2 + lab - b 2 ) . 


Or, 


«/,(a + ) 


- n 2 


- 


(6 - aY 


= a 2 - 


h 2 


This solution is a poor approximation to the exact solution u(a ) = a 2 . However, observe that at the Radau 
point (away from the right boundary b ), namely, t = a + h/ 3, the linear function u h takes on the value 

ui x (a + h/ 3) -(a + h/ 3) 2 , 

which is again exact. See Figure 4.1. The above observation is consistent with the next assertion. 


4,1 Equivalence Between DG and Radau IIA Methods 

If the right Radau quadrature is employed, then the DG method via (4.5) is identical to the collocation 
method using 5 right Radau points as collocation points, i.e., the Radau IIA method. (As a consequence, 
due to the equivalence result in the previous section, with such a quadrature, the DG, CG, and collocation 
methods become identical.) 

The key idea for the proof is the following. Since u h has a discontinuity or a jump across t„, the value 
u h '(t n ) is not well defined (it includes the derivative of a jump, which involves the Dirac delta function). 
Viewing this from another angle, the derivative u h '(t) on /„ incoiporates no effect from the known 
boundary value w/,(f„ - ) . As a result, u' for the ODE will be evaluated not by u h ' but by U where U, a 
function to be constructed, is continuous across all t n , and U (t n ) = Uj,(lY) ■ The continuous function U 
will be obtained by adding a correction to u h . 




Figure 4.1. — The piecewise linear (s = 2) DG solution for u'(t) = 2f with initial condition u(0) = 0 and h - 1. (a) After 
one step, the solution Uh is exact at the Radau points t = 1 and t = 1/3; (b) the solution after three steps (note the 
different scales). 
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To define U, we need the following correction function g, which serves the purpose of eliminating the 
trial function v (collocation methods do not use trial functions). This correction function, introduced for 
conservation laws in (Huynh 2006), is where the current approach differs from those in the literature. It 
leads to a DG formulation using the differential instead of integral equations, and results in a simplified 
version of the DG method. 

To deal with the correction term on the left hand side of (4.5), in view of the expression 

j* f ' 

Uh ( t)v(t ) dt , we ask the following question: can we define a polynomial g on /„ which possesses the 

t n 

property that for all v in P ( 2 ', 


J " +1 g'(t)v(t ) dt = - v(t„) . (4.7) 

Such a function will help eliminate the correction term and the test function. The answer to the above 
question is positive: applying integration by parts to the left hand side of (4.7), we have 

[ " +l g'(t)v(t) dt = g(t n+ i)v(t n+l ) - g(t n )v(t n ) - [ " +1 g{t)v'(t) dt. (4.8) 

J t n " t n 

Thus, for (4.7) to hold, it suffices that 

gifn) = 1, g(t„+ 1 ) = 0, (4.9) 


and 


[ , +1 g(t)v'(t) dt - 0 for all v in P ? 1 . (4.10) 

" tn 


Since v is in P s ', V is in P s 2 ; in addition, as v spans P s ', v' spans P 2 . Consequently, the above is 
equivalent to g being orthogonal to P : 


[ " +1 g(t)w(t) dt - 0 for all w in P' 5 ' 2 . (4.11) 


That is, 


f tn+1 g(t)(t-t n )J dt = 0 for j - 0,1,..., 5 - 2 . (4.12) 

" tn 

In other words, g approximates the function 0 in a least square sense. The requirement (4.1 1) or (4.12) 
provides s — 1 conditions to define g. Together with the two conditions in (4.9), we have a total of .v + 1 
conditions. Therefore, we require g to be of degree 5. 

On the interval /„, the right Radau polynomial of degree s satisfies conditions (4.1 1) and (4.9) (see 
Fig. 4.2(a)). Let g be this right Radau polynomial. Then g can also defined by g (t„) =1 and g vanishes at 
the s right Radau points on /„. (For the relation between Legendre and Radau polynomials, see the 
Appendix.) Next, set 


U — M/; + [W/;(6; ) Uj % (t n )]g. 


(4.13) 
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Figure 4.2. — The DG method with s = 3. (a) The correction function g is the Radau polynomial of degree s on 
/ = [0,1]. It is defined by: g(0) =1, g(1) = 0, , and the projection of g onto P s-2 is the zero function, (b) The 
polynomial u h is of degree s - 1 and is generally discontinuous across t n (thin curve); the polynomial U is of degree 
s and is continuous across f„ (thick curve); U is defined by: U(t n ) = u h (t n ) , U(t n+l ) = u h (t n ~ i ) , and the projections 
of U and u h onto P s ~ 2 are the same. 


Then, U is of degree 5 (whereas, u h is of degree 5 - 1) as shown in Figure 4.2(b). Using (4.9), we obtain 


U(t n ) w/j (t n ) 


(4.14a) 


and 


u 0 n+1 ) = u h0 n+l ). 


(4.14b) 


Since the above two equations hold for all n, replacing n + 1 by n in the latter, we have U(t„ ) = u h (t n ) • 
Therefore, by (4.14a), 


On ) — U(t n ) — U h On ) • 

That is, U is continuous across all t„. Next, by (4.13), 

U — U ] 7 + [UhOn ) — Uh On )]^ • 


As a consequence, using (4.7), 

[ " +l U'(t ) v(0 dt = - [u h 0n)-u h 0n)]v0n) + f " + ' U h ' 0)v(t) dt . 

J t n " t n 

Therefore, by (4.5), 

\ tn+X U'(t) v(0 dt = V n+l f(t,u h 0)) v(0 dt . (4.15) 

j t n j t n 

Note that there is no correction term in the above weak formulation. Also note that whereas U appears on 
the left hand side, u h appears on the right hand side. We need to replace «/, by U. To this end, let t n j, 
j = l,...,s, be the 5 right Radau points. Since g vanishes at these points, by (4.13), 

U0n,j) = Uh0 n j)- ( 4 -16) 

We now use the right Radau quadrature. If the weights of the quadrature are bj, then (4.15) implies 
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Yj b J H’ttnj) v(t nJ ) = Yjbj f(tn,j,U h {t n j)) v(t n j ) . 

j = 1 j = 1 

Let v be one of the basis functions L h which are the Lagrangian polynomials of degree .v - 1 defined by 
(2.8) for the right Radau points. Then 

bjU (tfi'i) — bj f (t n j,u/j(t n j)) . 


Therefore, by (4. 16), for i = 1 


U'(t nJ ) = f(t nJ ,U(t n j)). (4.17) 

The fact that U is of degree s, (4.14a), and the above shows that U is the solution by collocation method 
with the right Radau points as collocation points. This completes the proof. 

4,2 Remarks 

The following remarks are in order. 

The CG and DG solutions relate to each other as follows. Assume that the right Radau quadrature is 
employed, and let u n = u^it „ _ ) . For the CG method, the solution polynomial U CG is of degree 5 and 
determined by u n and the values «„|,...,«„ V at the .v right Radau points. For the DG method, the 
polynomial u h is of degree .v - 1 and determined by the values ..,u, us but not u„. See Figure 4.2(b). 

The polynomial U DG = U defined by (4.13) above is identical to U CG ■ Also note that with the right Radau 
quadrature, since c s = 1, we have t, us = t„+ \, and 

DG (f«+l) = u n,s = u h( t n + 0 = UcG( t n+l) ■ (4-18) 

Next, we discuss accuracy. Since u h is of degree s — 1, we expect u h (t) to have an error of 0(h ). On the 
other hand, U CG is of degree s ; therefore, its error on I„ is 0(h '). Note that u h and U CG take on the same 
values at t nG ,...,t, hS , namely, m„j,. respectively (see Fig. 4.2(b)). As a result, for the DG method, as 
values of u/,, the solutions at the interior right Radau points, namely, w„ i,. have errors 0(h '), one 

order higher than expected. As for u n s = u/ t (t~ + j) , it has an error of 0(h 2s ). For example, with 5 = 2, i.e., 
a piecewise linear the value u n 2 = u/,(t n ~ +l ) is accurate to 0(h 4 ), and the method is of order 3, 
whereas, u, h \ = u h (t n \) has an error of 0(h 3 ), i.e., it is as good as a parabolic approximation. 

The next remark concerns the function g defined by (4.7), i.e., for all v in P s_l , 

f " +1 g'(t)v(t) dt = - v(t n ) . 

J t n 

That is, for any v in P , the projection of v onto the line spanned by g' is the value —v(t n ). Using 
distribution theory, we can view g in the following manner. The above implies that on the subspace P'~', 
the function -g' is equivalent to the Dirac delta function at t = t„: 

-g'-S(G). 

Since the integral of the Dirac delta function is the unitstep (or Heaviside) function, 

f 0 if t < t n , 

-g + constant * H(t-t„) = < 

ll if t n <t<t n+ 1. 
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By choosing the constant to be 1, we obtain the result that g approximates the ‘step-down function’: 


[ 1 if t<t„, 

g « < 

{ 0 if t n <t<t n+l . 

That is, g(t „ ) = 1 , and g approximates the zero function on /„. 

4.3 A Generalization for DG Formulation 

The DG formulation above can be generalized. One way is to use, for each boundary t„, a quantity that 
blends the upwind value u/, (t n ~) and the downwind value u/, (t n + ) . Thus, for the left boundary 
evaluation in (4.3), the value to be employed is, with a a parameter between 0 and 1, 

u n = a u h {t n -) + (1 -a)u),it n + ) . (4.19) 

If a = 1, then u n = ui,(t ,,) ; a = 1/2, then u n = [w/, (hD + Uh(t n + )\/ 2; and, a = 0, then u n = Uf,(t n + ) . 
This approach to generalization was presented by Delfour, Hager, and Trochu (1981). For the right 
boundary t n + 1 , in addition to W/,(t n +j), we also need m/j(*„+j) . The use of such a value results in a 
complicated method. 

Here, we introduce a generalization that does not employ Uf,(t *j) by formulating the CG method in 
the form of DG. Suppose u„ is given, and u n + 1 is to be calculated. Consider the CG method using, for the 
moment, the Gauss quadrature whose evaluation points are the Gauss points denoted by t n \,. . ,,t ns where 
t n < t, L \ < ... t, hS < t„ i . As shown in Section 3.0, the solution polynomial U = Ucg is identical to the 
collocation solution with collocation points t n ,i,...,t„^. That is, U is of degree s and interpolates u n and the 
s collocation values ..,»„, s . Consider the polynomial of degree s — 1 interpolating «„ lv (but not 
u n ) denoted by u h and shown in Figure 4.3(b). This polynomial plays the same role as u h of the DG 
method. In addition, on I = [0, 1 ], let g be a polynomial of degree .v that takes on the value 1 at x = 0 and 
vanishes at t„j,. . . i.e., g is the function <t > 0 defined in (3.9) and shown in Figure 4.3(a) (here, g is the 

Legendre polynomial on [0,1]). The polynomial U = U C g can be written as, with t = t n + x h, 

U it) = u h it) + [u n - u h (f„ + )] gUt -t n )l h) . (4.20) 

The solution we are seeking is u n+ \ = U(t„ i), 

U n + 1 = U(t n+ 1 ) = Uit n ~) = u h (/„-)+ [u n -Uhitn)] g(l). (4.21) 

That is, u n +\ is obtained by adding to the value «/, it ri ~ l ) a correction term which depends only on the 
jump u n - Uf, (t„ + ) at the left boundary. The question is: if U is the CG solution, then what is the 
corresponding DG formulation for «/,? Using the test space P in the CG formulation, for any v in P s_1 , 

iU',v)=if,v). (4.22) 

Applying integration by parts to the left hand side above, 

r t fj+\ 

iu',v) = w„+iv(G +1 )-M„ v(G +1 )- U(t)v'it)dt. (4.23) 

J t n 


Next, by (4.20), 
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[ " +1 U(t)v'(t)dt = \" +X {u h (t)+ \u n -u h (t+)\ g((t-t n )l h)} v\t)dt . 

J t n " t n 

Using the fact that g is orthogonal to P 2 (in fact, g is orthogonal to P s ~'), 

[ ” +1 U(t) v'(t)dt — [ " +1 Ufo (t) v'(t)dt . 

* t n " t n 

where, again, v is in P‘ s_1 Thus, (4.23) implies 

C ^n+l 

(U',v) = u n+ i v(t„+i)-u n v(t n )-\ u h (t)v'(t)dt . (4.24) 

Applying integration by parts on the last term above, we obtain 

I u h (t)v'(t)dt = u h (t n ~ x )v(t n+ i)-u h (t n )v{t n )-\ u h (t)v(t)dt . 

J t n " t n 

The above two equations imply 

(U',v) = K+i-w/^+^Mtw+i) - [u„-u h (t n )\v{t n ) + u h (t)v(t)dt. (4.25) 
Next, by (4.21), 

“n+1 ~ u h 0„+i) = [tl n -Uh(fn)\ gC 1 ) • 

Substitute the above into (4.25), 

(U',v) = \u n -u h (t n )] g(l)v(t„ +1 ) - [u n -u h {t n )]v{t n ) + u h (t)v(t)dt. 

" t-n 

Or 


(U',v) = [u n -u h {tn)\ [g(l)v(t„+i) -v{t n )] + [ " +1 U h '(t) v(t)dt. (4.26) 

~ t-n 

Therefore, by (4.22), 

[tin -Uhitn)] [g(l)v(^+l) -v(t„)] + (u h ',v) = (/,v). (4.27) 

Note that the above reduces to the DG formulation (4.5) if g(l) = 0, i.e., if x, =1. 

This completes the generalization of the DG method. Also note that such a DG formulation is more 
involved than the CG formulation. 
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(a) 



Figure 4.3. — A generalization DG method for s = 3. (a) Here, the quadrature points are the Gauss points. The 
function g on / = [0,1] is of degree s and defined by the conditions that g(0) = 1, and g vanishes at the quadrature 
points; for this case, g is the Legendre polynomial of degree s. (b) The polynomial Uh is of degree s - 1 and 
interpolates the values u n ,i,...,u ni s at the s quadrature points. 


5.0 Conclusions and Discussion 

In summary, we studied the numerical solutions for ordinary differential equations using the 
collocation, continuous Galerkin (CG), and discontinuous Galerkin (DG) methods. It was shown that if a 
quadrature formula using s evaluation points is employed, then the CG method is identical to the 
collocation method using quadrature points as collocation points. Furthermore, if the quadrature formula 
is the right Radau one, then the DG and CG methods also become identical, and they reduce to the Radau 
IIA collocation method. We also present a generalization of DG that yields a method identical to CG and 
collocation with arbitrary collocation points. As a result of these findings, both the CG and DG methods 
can be formulated using the differential instead of integral form. 

In addition to clarifying the relations among these methods, the equivalence results can be employed 
for the high-order time discretization of conservation laws. It can also be applied to simplify the time 
discretization of the space-time DG methods. 
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Appendix A. — Radau Polynomials 

Since the Radau polynomials are not widely known, we derive them below. The Radau polynomial of 
degree s, which is determined by conditions (4.9) and (4.1 1) (or (4.12)), is defined here by using the 
Legendre polynomials. The advantage of such a definition is that it clarifies the relation between these 
polynomials as well as the orthogonality properties (to the various space of polynomials). Instead of the 
interval [0,1], to be consistent with the standard Legendre polynomials, we employ the interval [-1,1]. If 
g is a polynomial on [-1,1], then the corresponding polynomial on [0,1] can easily be obtained by, for q 
on [0,1], G(p) = g(2r\ - 1). 

We now focus on/= [-1,1]. For any integer m> 0, let P„, be the space of polynomials of degree m or 
less. Then P,„ is a vector space of dimension m + 1. A polynomial v is orthogonal to P,„ if, for each l, 
0<l<m, 


= o. 

Clearly, the criterion of being orthogonal to P„, provides m + 1 conditions (or equations). 

For k = 0,1,2,..., let the Legendre polynomial P k be defined on 1= [-1,1] as the unique polynomial of 
degree k satisfying the following k +1 conditions: it is orthogonal to P/,_| and P k (1) = 1. The Legendre 
polynomials are given by a recurrence formula (e.g., Hildebrand 1987): 

Po = 1, I\=t» 


and, for k> 2, 


Pk&) = lk T 1 1© - V P k~ 2(5)- (A- 1 ) 

k k 

The first few Legendre polynomials are plotted in Figure A. 1(a). Useful properties of the Legendre 
polynomials are listed below. If k > m, then P k is orthogonal to P,„. Next, P k is an even function (involving 
only even powers of Q for even k, and an odd function for odd k. For all k, the values at the boundaries 
are 


P k (-\) = (-\) k , (A.2a) 

Pk( 1) = 1. (A.2b) 


The derivative values at the end points are 


Pk = (-l) k ~ l k(k + l)/2 , (A. 3 a) 

Pk {\) = k(k + 1)/2 . (A.3b) 

In addition, for k±l, (P k J J i) = 0. Finally, the zeros of P k are the k Gauss points on [-1 , 1 ]. 

The right Radau polynomial of degree k (k > 1 ) is defined by 

RR,k = tj^(Pk-Pk-l). (A.4) 
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The letter R stands for ‘Radau’ and the subscript R for ‘right’. The factor (-1 / is nonstandard and is 
needed so that (A.6a) below holds. The first few Radau polynomials are plotted in Figure A. 1(b). The 
above definition implies that R R k is orthogonal to P/,_ 2 - In addition, by (A.2), 


%(-!)=!■ (A. 5 a) 

R R ,kV)= 0 (A.5b) 


It is important to note that R Rk , which is of degree k, is defined by the above two conditions and the k - 1 
conditions that it is orthogonal to P /c _ 2 . This definition of the Radau polynomial shows that it 
approximates the zero function in the sense of least squares. At the two boundaries, by using (A.3), 

R R , k '(-l) = -k 2 /2, (A. 6a) 

R R , k '(l) = (-l) k - l k/2. (A.6b) 


The zeros of the Radau polynomial R R j are the bright Radau points. 

For later use, the Lobatto polynomial of degree k (k > 1 ) is defined by 

Lo k = p k ~Pk-2 ■ 


(A.7) 


The zeros of the Lobatto polynomial of degree k are the k Lobatto points; they include the two boundaries 




- 0.5 


0 0.5 1 -1 - 0.5 0 

(a) (b) 

Figure A.1 . — (a) Legendre polynomials and (b) right Radau polynomials. 


0.5 
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Appendix B. — Examples of Collocation Methods 

In the following examples of collocation methods, we make several uncommon observations, which, 
hopefully, will shed some light on the matter. First, we write the collocation method in IRK form. Then 
we calculate the solution for the ODE u' = A .u by assuming that 

u n = 1 and h = 1. (B.l) 

The assumption h= 1 means, loosely put, h is absorbed into z where z = )Ji = h and, for order of accuracy 
calculation, halving the step size takes the form of halving z. As for the stability function, 

R(z ) = u n+ 1 . (B.2) 

The case of 2 Lobatto points corresponds to a Butcher array shown in Table B.l (a). The resulting 
collocation method reduces to the Trapezoidal Rule: 

Un+l = Un + 2 'h(fn fn + 1 ) • (B.3) 

The stability function and its error are, respectively, 

R(z) = (l + iz)/(l-iz),and E(z) = e-- - R(z) = -^-z 3 + 0(z 4 ). (B.4) 

Thus, the method using 2 Lobatto points is of order 2. It is also A-stable as can be seen by Figure B. 1(a). 


TABLE B.l.— BUTCHER ARRAYS FOR 2-POINT COLLOCATION METHODS. 


0 

0 

0 

1/2-V3/6 

1/4 

1/4 - V3/6 

1 

1/2 

1/2 

1/2 + V3/6 

1/4 + V3/6 

1/4 


1/2 

1/2 


1/2 

1/2 


(a) Lobatto (b) Gauss 


For the case of two Gauss points, the Butcher array is shown in Table B. 1(b). Using assumption (B.l), 


_ 1-W3/6 

U "’ 1 ~ l-A/2 + A 2 /12 


and u „2 


1 + W3/6 
l-A/2 + A 2 /12 ' 


The solution at t„+\ and the stability function are, respectively, 


Ufl + 1 


l + A/2 + A 2 /12 , . 1 + z/2 + z 2 /12 

and R(z) = 

l-A/2 + A 2 /12 1-z/2 + z 2 /12 


(B.5) 


(B.6) 


At the collocation points, u nJ and m„ 2 respectively approximate e z c ' and e zc ' 2 ; the errors are, 

e n \ = e ZCl - u n \ * 0.008019z 3 + 0.005167z 4 + 0.002008z 5 + G(z 6 ), and 
e n 2 = e ZC2 — u n 2 ~ - 0.008019z 3 - 0.005167z 4 + 0.002008z 5 + 0(z 6 ). 
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(a) 2 Lobatto points (b) 2 Gauss points 

Figure B.1 . — Plot of |R(z)| for 2-point collocation methods where z = x + iy varies on the complex plane. On the region 
{z; |R(z)| > 1}, the plot is cut off (flat part); the region z; |R(z)| < 1} is the stability domain S. The two methods are A- 
stable: (a) 2 Lobatto points and (b) 2 Gauss points. 


As expected, these errors are 0(z 3 ). The cancellation of errors yields a solution u„+\ with an error of 0(z 5 ): 

E(z) = e z - R(z) = z 5 / 720 + 0(z 6 ) = 0.001389z 5 + 0(z 6 ) . (B.8) 

Thus, the (2 -stage) collocation method using 2 Gauss points is of order 4. For comparison, the leading term 
of the error of the 4-stage explicit RK method is z s /5! = z 5 /120 = 0.008333z 5 , which is six times larger. 
Concerning stability, the method using two Gauss points is A-stable as can be seen in Figure B.l(b). 

Note that with two points (or two evaluations), the Gauss quadrature is exact for a cubic (its degree of 
precision is 3). Next, if the exact u' is a cubic, then the exact u is a quartic and, in this case, u n+ \ is exact. 
Since w„ +1 is exact for the case of a quartic u, for the general case, u n +\ has an error of 0(h 5 ), a fact 
consistent with (B.8). 

For the case of 2 left Radau points, the Butcher array is shown in Table B.2(a). Using assumption 
(B.l), 


u n \ = 1 and u n 2 = (1 + A, / 3)/(l — A, / 3) . 


(B.9) 


The solution at f„ +1 and the stability function are, respectively, 


Un + 1 


1 + 2A./3 + ^ 2 /6 
1 - A/3 


and R(z) 


l + 2z/3 + z 2 /6 
l-z/3 


At c 2 = 2/3, u „_2 approximates e“ 3 with an error of 


e n,2 = ^ 2z/ 3- M„ 2 =-(2/81)z 3 + 0(z 4 ) . 
The solution u„ +1 approximates e : with an error of 

E(z) = e z - R(z) = -z 4 /72 + 0(z 5 ) . 


(B-10) 


(B.l 1) 


(B.12) 


That is, u fl 2 has an error of G(z’), whereas u n+u 0(z 4 ). Thus, the collocation method using 2 left Radau 
points is of order 3. Concerning stability, this method is not A-stable as can be seen in Figure B.2(a). 
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(a) 2 left Radau points (b) 2 right Radau points 

Figure B.2. — Plot of |R(z)| for 2-point collocation methods where z = x + iy varies on the complex plane. On the region 
{z; |R(z)| > 1}, the plot is cut off; the region {z; |R(z)| < 1} is the stability domain S. (a) The left Radau method is not 
A-stable. (b) The right Radau method is L-stable, thus, also A-stable. 


TABLE B.2.— BUTCHER ARRAYS FOR 2-POINT COLLOCATION METHODS 


0 

0 

0 

1/3 

5/12 

-1/12 

2/3 

1/3 

1/3 

1 

3/4 

1/4 


1/4 

3/4 


3/4 

1/4 


(a) 2 left Radau points (b) 2 right Radau points 


For the case of 2 right Radau points, the Butcher array is shown in Table 2.2(b). Using assumption 
(B-l), 


^n, 1 


l-X/3 

1 - 2A / 3 + A 2 / 6 


and 


Mn, 2 


l + X/3 

1 - 2A / 3 + A 2 / 6 ‘ 


(B.13) 


Since c 2 = 1, we have u, u2 = u n+h and 


R(z) 


l + z/3 

l-2z/3 + z 2 / 6 


At ci = 1/3, u„_ i approximates e n with an error of 


e n i = e zl ^-u n i = (2/81)z 3 + 0(z 4 ) 


The solution w„+i approximates e : with an error of 


(B.14) 


(B.15) 


E(z) = e z -R{z) = e z -u n+ \ = z 4 /72 + 0(z 5 ) . (B.16) 

That is, «„,! has an error of (T(z’), whereas, u n+ \, 0(z 4 ). Thus, the (2 -stage) collocation method using 2 
right Radau points is of order 3. For comparison, the error of the standard 3-stage explicit Runge-Kutta 
(RK) method is A 4 /(4!) = A 4 /24, which is three times larger. Concerning stability, this Radau method is In- 
stable as can be seen in Figure B.2(b). Collocation methods using right Radau points are called Radau IIA 
methods (see, e.g., Lambert 1991 or Hairer and Wanner 1991). 

Note that the errors for the left Radau case in (B.l 1) and (B.12) and those for the right Radau case in 
(B.15) and (B.16) are of opposite sign and same magnitude. Also note that the curve {z; \R(z)\ = 1 } for the 
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left Radau method in Figure B.2(a) and that for the right Radau method in Figure B.2(b) are reflections of 
each other about the origin of the complex plane. Indeed, the equations for these curves are, respectively, 
by (B. 10) and (B. 14), 


|l + jz + |z 2 | = |l-±z| and 1 -fz + ±z 2 | = 1 1 + jz | . (B.17) 

Replacing z by —z for the equation on the left, we obtain that on the right, and the observation follows. 

The above observations on the errors as well as the symmetry of the curves {z; \R(z)\ = 1 } between the left 
and right Radau methods hold for the general case of s collocation points as well. 

Finally, loosely put, the more biased toward t „+ 1 are the collocation points (e.g., right Radau points), 
the more stable is the method. 
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